% !TEX root = FDS_Technical_Reference_Guide.tex

\typeout{new file: Radiation_Chapter.tex}

\chapter{Thermal Radiation}
\label{chapter:radiation}

Gas phase thermal conduction and radiation are represented by the divergence of the heat flux vector in the energy equation, $\nabla \cdot \dbq''$. This chapter describes the equations that govern the radiative component, $\dbq''_{\rm r}$.


\section{Radiation Transport Equation}
\label{RTEsection}

The Radiative Transport Equation (RTE) for an absorbing, emitting, and scattering medium is \cite{Viskanta:1987}
\begin{align}
\bs \cdot \nabla I_{\la}(\bx,\bs) &=
\underbrace{ - \kappa(\bx,\la)   \; I_\la(\bx,\bs) }_{\textrm{Energy loss by absorption}} -
\underbrace{\sigma_{\rm s}(\bx,\la) \; I_\la(\bx,\bs) }_{\textrm{Energy loss by scattering}} +  \notag \\[0.25in]
& \underbrace{   B(\bx,\la) }_{\textrm{Emission source term}} + \quad
\underbrace{   \frac{\sigma_{\rm s}(\bx,\la)}{4\pi}
\int_{4\pi}\Phi(\bs',\bs) \; I_{\la}(\bx,\bs') \; d\bs'
 }_{\textrm{In-scattering term}}
\label{RTEbasic}
\end{align}
where $I_{\la}(\bx,\bs)$ is the radiation intensity at wavelength, $\la$; $\bs$ is the direction vector of the intensity; and $\kappa(\bx,\la)$ and $\sigma_{\rm s}(\bx,\la)$ are the local absorption and scattering coefficients, respectively. $B(\bx,\la)$ is the emission source term, describing how much heat is emitted by the local mixture of gas, soot and droplets/particles. The integral on the right hand side describes the in-scattering from other directions. The in-scattering and scattering terms are detailed in Section~\ref{droplet-radiation}.

In practical simulations, the spectral dependence of the RTE cannot be resolved accurately. Instead, the radiation spectrum is divided into a relatively small number of bands and a separate RTE is derived for each band. For instance, the band specific RTE for a non-scattering gas is
\be
   \bs \cdot \nabla I_n(\bx,\bs) = B_n(\bx) - \kappa_n(\bx) \, I_n(\bx,\bs),\;\; n = 1...N  \label{bandRTE}
\ee
where $I_n$ is the intensity integrated over the band $n$, and $\kappa_n$ is the appropriate mean absorption coefficient for the band. When the intensities corresponding to the bands are known, the total intensity is calculated by summing over all the bands
\be
   I(\bx,\bs) = \sum_{n=1}^N I_n(\bx,\bs)
\ee

\subsection{Radiation Source Term}

The emission source term for radiation band $n$ is
\be
   B_n(\bx) = \kappa_n(\bx) \, I_{{\rm b},n}(\bx)
\ee
where $I_{{\rm b},n}$ is the fraction of the blackbody radiation at temperature $T(\bx)$:
\be
   I_{{\rm b},n}(\bx) = F_n(\la_{\rm min},\la_{\rm max}) \, \sigma \, T(\bx)^4/\pi
\ee
and $\sigma$ is the Stefan-Boltzmann constant. The calculation of factors $F_n$ is explained in Ref.~\cite{Siegel:1}. The measurement of the absorption coefficients, $\kappa_n$, is discussed in Appendix~\ref{absorption_coefficients}.

Even with a reasonably small number of bands, solving multiple RTEs is very time consuming. Fortunately, in most large-scale fire scenarios soot is the most important combustion product controlling the thermal radiation from the fire and hot smoke. As the radiation spectrum of soot is continuous, it is possible to assume that the gas behaves as a gray medium.  The spectral dependence is then lumped into one absorption coefficient ($N=1$) and the source term is given by the blackbody radiation intensity \cite{Ludwig:NASA}
\be
   I_{\rm b}(\bx) = \frac{\sigma \, T(\bx)^4}{\pi} \label{emission_source_term}
\ee
This is the default mode of FDS. For optically thin flames, however, where the yield of soot is small compared to the yields of $\rm CO_2$ and water vapor, the gray gas assumption can lead to an over-prediction of the emitted radiation. From a series of numerical experiments using methane as the fuel, it has been found that six bands ($N=6$) provide an accurate representation of the most important radiation bands of the fuel, $\rm CO_2$, and water vapor~\cite{Hostikka:3}. Table~\ref{banditos} through Table~\ref{band_MMA} list the band limits for various fuel species. The location of the bands have been adjusted to accommodate most of the features of the fuels spectra. If the absorption of the fuel is known to be important, separate bands can be reserved for fuel, increasing the total number of bands, $N$. The number of additional bands depends on the fuel, as discussed in Appendix~\ref{absorption_coefficients}.

\begin{table}[p]
\caption{Limits of the spectral bands for methane (CH$_4$).}
\label{banditos}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 3400}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 2400}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 2174}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 1000}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                     & Soot           & CO$_2$       & CH$_4$ & CO$_2$ & H$_2$O,CH$_4$& H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot   & Soot   & Soot         & CO$_2$ \\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in} 1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 2.94}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 4.17}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 4.70}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in} 10.0}c@{200}}{ } \\
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Limits of the spectral bands for ethane (C$_2$H$_6$).}
\label{band_Ethane}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3350}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2550}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1650}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1090}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_2$H$_6$ & CO$_2$ & C$_2$H$_6$ & H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot  & CO, H$_2$O, Soot & H$_2$O, Soot & CO$_2$, C$_2$H$_6$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.99}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.92}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}6.06}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}9.17}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE ETHYLENE: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for ethylene (C$_2$H$_4$).}
\label{band_Ethylene}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3375}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1650}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}780}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_2$H$_4$ & CO$_2$ & C$_2$H$_4$ & H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot & CO, H$_2$O, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.96}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.57}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}6.06}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}12.82}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE PROPYLENE: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for propylene (C$_3$H$_6$).}
\label{band_Propylene}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3250}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2600}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1950}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1175}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_3$H$_6$ & CO$_2$ & C$_3$H$_6$ & C$_3$H$_6$, H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot & CO, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.08}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.85}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}5.13}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}8.51}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE PROPANE: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for propane (C$_3$H$_8$).}
\label{band_Propane}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3350}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2550}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1650}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1175}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_3$H$_8$ & CO$_2$ & C$_3$H$_8$ & H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot  & CO, H$_2$O, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.99}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.92}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}6.06}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}8.51}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE HEPTANE: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for heptane (C$_7$H$_{16}$).}
\label{band_Heptane}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3250}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2550}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1775}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1100}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_7$H$_{16}$ & CO$_2$ & C$_7$H$_{16}$ & H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot  & CO, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.08}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.92}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}5.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}9.09}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE TOLUENE: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for toluene (C$_7$H$_8$).}
\label{band_Toluene}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3200}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2550}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2050}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1200}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CO$_2$ & C$_7$H$_8$ & CO$_2$  & C$_7$H$_8$ & C$_7$H$_8$, H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & H$_2$O, Soot & Soot  & CO, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.12}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.92}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}4.88}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}8.33}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE METHANOL: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for methanol (CH$_3$OH).}
\label{band_Methanol}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3825}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3200}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2600}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1675}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1125}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot & CH$_3$OH & CH$_3$OH & CO$_2$ & CH$_3$OH & CH$_3$OH, H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species} & CO$_2$, H$_2$O & CO$_2$, Soot & Soot  & CO, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.61}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.12}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.85}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}5.97}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}8.89}c@{200}}{ } \\

\end{tabular}
\end{center}
\end{table}

%%% TABLE MMA: CHECKED ON 04/09/13

\begin{table}[p]
\caption{Limits of the spectral bands for methyl methacrylate (MMA).}
\label{band_MMA}
\begin{center}
\begin{tabular}{|*{7}{c|}}
\multicolumn{1}{c}{$\omega$ (1/cm)}
             & \multicolumn{1}{@{\hspace{-.2in}10000}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3800}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3250}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2650}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2050}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}1250}c@{50}}{ } \\
\hline
\hspace{0.2in} \underline{6 Band Model} \hspace{0.2in} & 1  & 2  & 3 & 4  & 5 & 6  \\ \cline{2-7}
                                      & Soot           & CO$_2$ & MMA & CO$_2$ & MMA & MMA, H$_2$O \\
\raisebox{1.5ex}[0pt]{Major Species}  & CO$_2$, H$_2$O & H$_2$O, Soot & Soot   & CO, Soot & H$_2$O, Soot & CO$_2$\\ \hline
\multicolumn{1}{c}{$\la$ ($\mu$m)}
             & \multicolumn{1}{@{\hspace{-.2in}1.00}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}2.63}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.08}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}3.77}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}4.88}c@{}}{ }
             & \multicolumn{1}{@{\hspace{-.2in}8.00}c@{200}}{ } \\

\end{tabular}
\end{center}
\vspace{2in}
\end{table}

\subsection{Radiation Contribution to Energy Equation}

The radiant heat flux vector $\dbq_{\rm r}''$ is defined
\be \dbq_{\rm r}''(\bx) = \int_{4\pi} \; \bs' \, I(\bx,\bs') \; \d \bs'   \ee
The gas phase contribution to the radiative loss term in the energy equation is (under the gray gas assumption)
\be -\nabla\!\cdot \dbq_{\rm r}''(\bx)(\mbox{gas}) =
    \kappa(\bx) \, \left[ U(\bx) - 4 \pi \, I_{\rm b}(\bx) \right]  \quad ; \quad
    U(\bx) = \int_{4\pi} \, I(\bx,\bs') \, \d \bs'  \label{net_emission}
\ee
For N bands, the gas phase contribution to the radiative loss term in the energy equation is
\be
-\nabla\!\cdot \dbq_{\rm r}''(\bx)(\mbox{gas}) = \sum_{n=1}^N  \kappa(\bx) \, U_n(\bx) \, -4 \pi B_n(\bx) \quad
; \quad U_n(\bx) = \int_{4\pi} \, I_n(\bx,\bs') \, \d \bs'  \label{net_emission_Nbands}
\ee

In words, the net radiant energy gained by a grid cell is the
difference between that which is absorbed and that which is emitted.

\subsection{Correction of the Emission Source Term}

In calculations of limited spatial resolution, the source term, $I_{\rm b}$, defined in Eq.~(\ref{emission_source_term}) requires special treatment in the flaming region of the fire. Typical FDS calculations use grid cells that are tens of centimeters in size, and consequently the computed temperatures constitute a bulk average for a given grid cell and are considerably lower than the maximum temperature in a diffusion flame. Because of its fourth-power dependence on the temperature, the source term must be modeled in those grid cells where combustion occurs. Elsewhere, the computed temperature is used directly to compute the source term. It is assumed that this ``flaming region'' is where the local, nominal radiative loss is greater than a specified lower bound, $\chi_{\rm r} \dq'''>10$~kW/m$^3$. In this region, the global radiative fraction model is used. The emission source term is multiplied by a corrective factor, $C$:
\be I_{\rm b,f}(\bx) = C \, \frac{\sigma \, T(\bx)^4}{\pi}   \quad ; \quad
    C = \min \left( 100 \; , \; \max \left[ 1 \; , \; \frac{\sum_{\dq'''_{ijk}>0} \left( \chi_{\rm r} \, \dq'''_{ijk} + \kappa_{ijk} \, U_{ijk} \right) \, \d V}{\sum_{\dq'''_{ijk}>0} \left( 4 \, \kappa_{ijk} \, \sigma \, T^4_{ijk} \right) \, \d V} \right] \right) \label{corrected_emission_source_term}
\ee
When the source term defined in Eq.~(\ref{corrected_emission_source_term}) is substituted into Eq.~(\ref{net_emission}), the net radiative emission from the flaming region becomes the desired fraction of the total heat release rate. Note that this correction factor is bounded below by 1 and above by 100. The upper bound is arbitrary, meant to prevent spurious behavior at the start of a simulation.

The radiative fraction, $\chi_{\rm r}$, is a useful quantity in fire science. It is the nominal fraction of the combustion energy that is emitted as thermal radiation. For most combustibles, $\chi_{\rm r}$ is between 0.3 and 0.4~\cite{Beyler2:SFPE}. However, in Eq.~(\ref{corrected_emission_source_term}), $\chi_{\rm r}$ is interpreted as the fraction of energy radiated from the flaming region.  For a small fire with a base diameter less than approximately 1~m, the local $\chi_{\rm r}$ is approximately equal to its global counterpart. However, as the fire increases in size, the global value will typically decrease due to a net re-absorption of the thermal radiation by the increasing smoke mantle~\cite{Takahashi:1}.

To account for the possibility that different fuels may have different values of radiative fraction, $\chi_{\rm r}$ is defined on a per-reaction basis, similar to the approach discussed by Gupta et al.~\cite{Gupta:2015}. If multiple $\chi_{\rm r}$ values are given, FDS will generate a local $\chi_{\rm r}$ by weighting the reaction specific $\chi_{\rm r}$ values by the local reaction rates.


\section{Numerical Method}
\label{radnumericalmethodsection}

This section describes how $\nabla\!\cdot \dbq_{\rm r}''$ (the radiative loss
term) is computed at all gas-phase cells, and how the
the radiative heat flux $\dq_{\rm r}''$ is computed at solid boundaries.

\subsection{Angular Discretization}
\label{radiation-discre}

To obtain the discretized form of the RTE, the unit sphere is divided into a finite number of solid angles.
The coordinate system used to discretize the solid angle is
shown in Figure~\ref{Angular}.
\begin{figure}[ht]
\begin{center}
\includegraphics[height=2in]{FIGURES/RadCoord}
\caption{Coordinate system of the angular discretization.}
\label{Angular}
\end{center}
\end{figure}
The discretization of the solid angle is done by dividing first
the polar angle, $\theta$, into $N_{\theta}$ bands, where
$N_{\theta}$ is an even integer.
Each $\theta$-band is then divided into
$N_{\phi}(\theta)$ parts in the azimuthal ($\phi$) direction.
$N_{\phi}(\theta)$ must be divisible by 4.
The numbers $N_{\theta}$ and $N_{\phi}(\theta)$ are chosen
to give the total number of angles $N_{\Omega}$ as close to
the value defined by the user as possible.
$N_{\Omega}$ is calculated as
\be
 N_{\Omega} = \sum_{i=1}^{N_{\theta}} N_{\phi}(\theta_i)
\ee
The distribution of the angles is based on empirical rules that try
to produce equal solid angles $\delta \bO^l = 4\pi/N_{\Omega}$. The
number of $\theta$-bands is
\be
 N_{\theta} = 1.17 \; N_{\Omega}^{1/2.26}
\ee
rounded to the nearest even integer. The number of $\phi$-angles
on each band is
\be
 N_{\phi}(\theta) = \max\left\{4,
        0.5\,N_{\Omega}\,\left[\cos(\theta^-)-\cos(\theta^+)\right]\right\}
\ee
rounded to the nearest integer that is divisible by 4.
$\theta^-$ and $\theta^+$ are
the lower and upper bounds of the $\theta$-band, respectively.
The discretization is symmetric with respect to the planes $x=0$, $y=0$, and
$z=0$. This symmetry has three important benefits:
First, it avoids the problems caused by the fact that the first-order
upwind scheme, used to calculate intensities on the cell boundaries,
is more diffusive in non-axial directions than axial.
Second, the treatment of the mirror boundaries becomes very simple, as
will be shown later. Third,
it avoids so called
``overhang'' situations, where $\bs\cdot {\bf i}$, $\bs\cdot {\bf j}$
or $\bs\cdot {\bf k}$ changes sign inside
the control angle. These ``overhangs'' would make the resulting system of
linear equations more complicated.

In the axially symmetric case these ``overhangs'' cannot be avoided, and a
special treatment, developed by Murthy and Mathur~\cite{Murthy}, is
applied. In these cases $N_{\phi}(\theta_i)$ is kept constant, and
the total number of angles is $N_{\Omega} = N_{\theta} \times
N_{\phi}$. In addition, the angle of the vertical slice of the cylinder is
chosen to be the same as $\delta\phi$.



\subsection{Spatial Discretization}


The grid used for the RTE solver is the same as for the hydrodynamic solver. The radiative transport equation~(\ref{bandRTE}) is solved using
techniques similar to those for convective transport in finite volume methods for fluid flow~\cite{Raithby}; thus, the name given to it is
the Finite Volume Method (FVM). More details of the model implementation are included in Ref.~\cite{Hostikka:2008}.

The thermal radiation spectrum is first divided into bands, as described in section~\ref{RTEsection}.  The procedure outlined below is
applied for each band of a wide band model, and thus the subscript $n$ has been removed for clarity.
The unit sphere is then discretized, as described in section~\ref{radiation-discre}.
Finally, the computational domain is divided into numerical grid as described in section~\ref{govequations}.
In each grid cell, a discretized equation is derived by integrating Eq.~(\ref{bandRTE}) over the volume of cell $ijk$ and the control
angle $\delta \Omega^l$, to obtain
\be
  \int_{\delta \Omega^l} \int_{V_{ijk}}
   \bs' \cdot \nabla I(\bx',\bs') \d\bx' \d\bs' =
   \int_{\delta \Omega^l} \int_{V_{ijk}} \kappa(\bx') \;
    \left[ I_{\rm b}(\bx') - I(\bx',\bs') \right] \d \bx' \d \bs'
\ee
The volume integral on the left hand side is replaced by a surface integral over the cell faces using the divergence theorem.
\be
  \int_{\delta \Omega^l} \int_{A_{ijk}}
   \left(\bs' \cdot \bn' \right) I(\bx',\bs') \d\bn' \d\bs' =
   \int_{\delta \Omega^l} \int_{V_{ijk}} \kappa(\bx') \;
    \left[ I_{\rm b}(\bx') - I(\bx',\bs') \right] \d \bx' \d \bs'
\ee
Assuming that the radiation intensity $I(\bx,\bs)$ is constant on each
of the cell faces, the surface integral can be approximated by a sum
over the cell faces.  Assuming further that $I(\bx,\bs)$ is constant within the volume $V_{ijk}$ and over the angle $\delta \bO^l$, and that
$\kappa(\bx')$ and $I_{\rm b}(\bx')$ are constants within the volume $V_{ijk}$, we obtain
\be  \sum_{m=1}^6 A_m \; I_m^l \;
      \int_{\Omega^l}(\bs' \cdot \bn_m) \d \bs'
   = \kappa_{ijk} \,
     \left[ I_{{\rm b},ijk} - I_{ijk}^l \right] \; V_{ijk} \,
     \delta \Omega^l   \label{RTEdiscrete2}
\ee
where
\begin{tabbing}
$I_{ijk}^l$ \hspace{1in}  \=  radiant intensity in direction $l$ \\
$I_m^l$                   \>  radiant intensity at cell face $m$ \\
$I_{b,ijk}$               \>  radiant blackbody Intensity in cell \\
$\delta \Omega^l$         \>  solid angle corresponding to direction $l$ \\
$V_{ijk}$                 \>  volume of cell $ijk$ \\
$A_m$                     \>  area of cell face $m$ \\
$\bn_m$                   \>  unit normal vector of the cell face $m$
\end{tabbing}
Note that while the intensity is assumed constant within
the angle $\delta \bO^l$, its direction covers the angle $\delta \bO^l$
exactly.The local incident radiation intensity is
\be
 U_{ijk} = \sum_{l=1}^{N_{\Omega}} I_{ijk}^l \delta\Omega^l
\ee

In Cartesian coordinates\footnote{In the axisymmetric case equation~(\ref{RTEdiscrete2}) becomes
a little bit more complicated, as the cell face normal vectors $\bn_m$ are not always constant. However, the computational efficiency can still be
retained.},
the normal vectors $\bn_m$ are the base vectors of the coordinate system and the integrals over the solid
angle do not depend on the physical coordinate, but the direction only. These integrals are denoted as
\be
D_m^l= \int_{\Omega^l}(\bs' \cdot \bn_m) \d \bs'
\label{cellfaceintegral}
\ee
and the discrete equation becomes
\be  \sum_{m=1}^6 A_m \; I_m^l \; D_m^l
   = \kappa_{ijk} \,
     \left[ I_{{\rm b},ijk} - I_{ijk}^l \right] \; V_{ijk} \,
     \delta \Omega^l   \label{RTEdiscrete3}
\ee
The cell face intensities, $I_m^l$ appearing on the left hand side of (\ref{RTEdiscrete3}) are calculated using a first order
upwind scheme. Consider, for example, a control angle having a direction vector $\bs$. If the radiation is traveling in the positive
$x$-direction, i.e., $\bs\cdot {\bf i} \geq 0$, the intensity on the upwind side, $I_{xu}^l$ is assumed to be
the intensity in the neighboring cell, $I_{i-1\,jk}^l$, and the intensity on the downwind side is the (unknown) intensity in the cell
itself $I_{ijk}^l$. The discrete RTE can now be written using the upwind intensities $I^l_{xu}$, $I^l_{yu}$ and $I^l_{zu}$ and $I_{ijk}^l$:
\begin{align}
   \notag A_x \; I_{xu}^l \; D_{xu}^l  & + A_x \; I_{ijk}^l \; D_{xd}^l +  \\
   \notag A_y \; I_{yu}^l \; D_{yu}^l  & + A_y \; I_{ijk}^l \; D_{yd}^l +  \\
   \notag A_z \; I_{zu}^l \; D_{zu}^l  & + A_z \; I_{ijk}^l \; D_{zd}^l  \\
   & = \kappa_{ijk}\,I_{{\rm b},ijk} V_{ijk} \,\delta \Omega^l
   - \kappa_{ijk}\,I_{ijk}^l       V_{ijk} \,\delta \Omega^l
\label{RTEdiscrete4}
\end{align}
where the $D$-terms on the LHS are integrals (\ref{cellfaceintegral}) evaluated on upwind and downwind sides of the cell. In the rectilinear mesh,
$D_{xu}^l = -D_{xd}^l$ and the equation can be simplified further. In addition, the integrals over the solid angle can be calculated analytically
\begin{align}
   D_x^l & = \int_{\Omega^l}(\bs^l\cdot {\bf i}) \, \d\Omega = \int_{\dph}\int_{\dth} (\bs^l\cdot{\bf i}) \sin\theta \; \d \theta \d\phi = \int_{\dph}\int_{\dth} \cos\phi\sin\theta\sin\theta\; \d\theta \d\phi \notag \\
         & = \frac{1}{2}\left( \sin\phi^+ - \sin\phi^- \right) \left[\Delta\theta - \left(\cos\theta^+\sin\theta^+ - \cos\theta^-\sin\theta^-\right)\right]  \\
   D_y^l & = \int_{\Omega^l}(\bs^l\cdot {\bf j}) \, \d\Omega  =  \int_{\dph}\int_{\dth} \sin\phi\sin\theta\sin\theta\; \d\theta \d\phi \notag \\
         & = \frac{1}{2}\left( \cos\phi^- - \cos\phi^+ \right) \left[\Delta\theta - \left(\cos\theta^+\sin\theta^+ - \cos\theta^-\sin\theta^-\right)\right]  \\
   D_z^l & = \int_{\Omega^l}(\bs^l\cdot {\bf k}) \d \Omega = \int_{\dph}\int_{\dth} \cos\theta\sin\theta\;  \d\theta \d\phi \notag \\
         & = \frac{1}{2}\Delta\phi \left[\left(\sin\theta^+\right)^2 - \left(\sin\theta^-\right)^2\right]
\\
   \delta \Omega^l & = \int_{\Omega^l} \d \Omega = \int_{\dph}\int_{\dth} \sin\theta \; \d\theta \d\phi
\end{align}
Equation (\ref{RTEdiscrete4}) for the unknown intensity $I_{ijk}^l$ is written in the form
\be
  a_{ijk}^l I_{ijk}^l =
  a_{x}^l   I_{xu}^l +
  a_{y}^l   I_{yu}^l +
  a_{z}^l   I_{zu}^l +   b_{ijk}^l \label{RTEdiscrete5}
\ee
where
\begin{align}
  a_{ijk}^l & = A_x |D_x^l| + A_y |D_y^l| +A_z |D_z^l| + \kappa_{ijk} \; V_{ijk} \delta\Omega^l \\
  \notag \\
  a_{x}^l & = A_x |D_x^l| \\
  a_{y}^l & = A_y |D_y^l| \\
  a_{z}^l & = A_z |D_z^l| \\
  b_{ijk}^l & = \kappa_{ijk} \; I_{{\rm b},ijk} \; V_{ijk} \; \delta\Omega^l
\end{align}
Here $\bf i$, $\bf j$ and $\bf k$ are the base vectors of the Cartesian coordinate system. $\theta^+$, $\theta^-$, $\phi^+$ and
$\phi^-$ are the upper and lower boundaries of the control angle in the polar and azimuthal directions, respectively, and $\Delta\theta =
\theta^+ - \theta^-$ and $\Delta\phi = \phi^+ - \phi^-$.

The solution method of (\ref{RTEdiscrete4}) is based on an explicit marching sequence~\cite{Kim}. The marching direction depends on the
propagation direction of the radiation intensity. As the marching is done in the ``downwind'' direction, the ``upwind'' intensities in all
three spatial directions are known, and the intensity $I_{ijk}^l$ can be solved directly from an algebraic equation. In the first cell to be solved,
all the upwind intensities are determined from solid or gas phase boundaries. In theory, iterations are needed if the reflections
or scattering are important, or if the scenario is optically very thick. Currently, no iterations are made.

\subsection{Boundary Conditions}

The boundary condition for the radiation intensity leaving
a gray diffuse wall is given as
\be I_{\rm w}(\bs) = \frac{\epsilon \, \sigma \, T_{\rm w}^4}{\pi} + \frac{1-\epsilon}{\pi}
 \int_{\bs'\cdot \bn_{\rm w} < 0} I_{\rm w}(\bs')\; |\bs'\cdot \bn_{\rm w} | \; \d\bs'
 \label{RTEbc} \ee
where $I_{\rm w}(\bs)$ is the intensity at the wall, $\epsilon$ is the
emissivity, and $T_{w}$ is the wall surface temperature.
In discretized form, the boundary condition on a solid wall is given as
\be I_{\rm w}^l = \frac{\epsilon \, \sigma \, T_{\rm w}^4}{\pi} + \frac{1-\epsilon}{\pi} \sum_{D_{\rm w}^{l'}<0} I_{\rm w}^{l'}\; |D_{\rm w}^{l'} |  \ee
where $D_{\rm w}^{l'}= \int_{\Omega^{l'}}(\bs\cdot \bn_{\rm w})\d\Omega$.
The constraint $D_{\rm w}^{l'}<0$ means that only the ``incoming'' directions
are taken into account when calculating the reflection.
The {\em net} radiative heat flux on the wall is
\be \dq_{\rm r}'' = \sum_{l=1}^{N_{\Omega}} I_{\rm w}^l \int_{\delta \Omega^l} (\bs' \cdot \bn_{\rm w}) \, \d\bs'
     = \sum_{l=1}^{N_{\Omega}} I_{\rm w}^l D_n^l \label{qrdef} \ee
where the coefficients $D_n^l$ are equal to $\pm D_x^l$, $\pm D_y^l$ or
$\pm D_z^l$, and can be calculated for each wall element at the start of the
calculation.

The open boundaries are treated as black walls, where the incoming intensity is
the blackbody intensity of the ambient temperature. On mirror
boundaries the intensities leaving the wall
are calculated from the incoming intensities using a
predefined connection matrix:
\be  I_{{\rm w},ijk}^l = I^{l'} \ee
Computationally intensive integration over all the incoming directions
is avoided by keeping the solid angle discretization symmetric on the $x$, $y$ and $z$ planes.
The connection matrix associates one incoming direction $l'$ to each mirrored direction on each wall cell.


\clearpage

\section{Absorption and Scattering by Small Droplets and Particles}
\label{droplet-radiation}

The attenuation of thermal radiation by liquid droplets and small particles is an important consideration, especially for water mist systems~\cite{Ravigururajan:1}.  Droplets and particles attenuate thermal radiation through a combination of scattering and absorption~\cite{Tuntomo:1}.  The radiation-spray interaction must therefore be solved for both the accurate prediction of the radiation field and for the particle energy balance. If the gas phase absorption and emission in Eq.~(\ref{RTEbasic}) are temporarily neglected for simplicity, the radiative transport equation becomes
\be
\bs \cdot \nabla I_{\la}(\bx,\bs) = -\Big[\kappa_{\rm p}(\bx,\la) + \sigma_{\rm p}(\bx,\la)\Big]
I_{\la}(\bx,\bs) +\kappa_{\rm p}(\bx,\la) \; I_{{\rm b,p}}(\bx,\la) +
\frac{\sigma_{\rm p}(\bx,\la)}{4\pi}
\int_{4\pi}\Phi(\bs,\bs') \; I_{\la}(\bx,\bs') \, \d\bs'
\label{RTEspray}
\ee
where $\kappa_{\rm p}$ is the particle absorption coefficient, $\sigma_{\rm p}$ is the particle scattering coefficient and $I_{{\rm b,p}}$ is the emission term of the particles. $\Phi(\bs',\bs)$ is a scattering phase function that gives the scattered intensity fraction from direction $\bs'$ to $\bs$.

\subsection{Absorption and Scattering Coefficients}

Radiation absorption and scattering by particles depends on their cross sectional areas and radiative material properties. For simplicity, we assume that the particles are spherical in shape, in which case the cross sectional area of a particle is $\pi r^2$, where $r$ is the particle radius. If the local number density distribution of particles at location $\bx$ is denoted by $n(r(\bx))$, the local absorption and scattering coefficients within a spray/particle cloud can be calculated from:
\begin{align}
\kappa_{\rm p}(\bx,\la) &= \int_0^\infty n(r(\bx)) \; Q_{\rm a}(r,\la) \; \pi r^2 \; \d r \\[0.1in] %N(\bx)\int_0^\infty f(r,d_m(\bx)) \; C_{\rm a}(r,\la) \; dr \\
\sigma_{\rm p}(\bx,\la) &= \int_0^\infty n(r(\bx)) \; Q_{\rm s}(r,\la) \; \pi r^2 \; \d r           %N(\bx)\int_0^\infty f(r,d_m(\bx)) \; C_{\rm s}(r,\la) \; dr
\end{align}
where $Q_{\rm a}$ and $Q_{\rm s}$ are the absorption and scattering efficiencies, respectively. For the computation of the spray/particle cloud radiative properties, spherical particles are assumed and the radiative properties of the individual particles are computed using Mie theory.

Based on Refs.~\cite{Collin} and \cite{Maruyama}, the real particle size distribution inside a grid cell is modeled as a mono-disperse suspension whose diameter corresponds to the Sauter mean ($d_{32}$) diameter of the poly-disperse spray. This assumption leads to a simplified expression of the radiative coefficients
\begin{align}
\kappa_{\rm p}(\bx,\la) &= A_{\rm p}(\bx) \; Q_{\rm a}(r_{32},\la) \\[0.1in]
\sigma_{\rm p}(\bx,\la) &= A_{\rm p}(\bx) \; Q_{\rm s}(r_{32},\la)
\end{align}
These expressions are functions of the total cross sectional area per unit volume of the droplets, $A_{\rm p}$, which is computed simply by summing the cross sectional areas of all the droplets within a cell and dividing by the cell volume. For practical reasons, a relaxation factor of 0.5 is used to smooth slightly the temporal variation of $A_{\rm p}$.

\subsection{Approximating the In-Scattering Integral}

An accurate computation of the in-scattering integral on the right hand side of Eq.~(\ref{RTEspray}) would be extremely time consuming and require a prohibitive amount of memory because the individual intensities in each location would have to be stored. Instead, a simplified form of Eq.~(\ref{RTEspray}) can be derived:
\be
\bs \cdot \nabla I_{\la}(\bx,\bs) =
-\left[\kappa_{\rm p}(\bx,\la) + \overline{\sigma}_{\rm p}(\bx,\la)\right]
I_{\la}(\bx,\bs)
+\kappa_{\rm p}(\bx,\la) \; I_{\rm b,p}(\bx,\la)
+\frac{\overline{\sigma}_{\rm p}(\bx,\la)}{4\pi}U(\bx,\la) \label{simple_particle_RTE}
\ee
where $U(\bx)$ is the total intensity integrated over the unit sphere and $\overline{\sigma}_{\rm p}$ is an effective scattering coefficient. The derivation of Eq.~(\ref{simple_particle_RTE}) is given in Appendix~\ref{radiation_derivations}. This simplified equation can be integrated over the spectrum to get the band specific RTE's. The procedure is exactly the same as that used for the gas phase RTE. After the band integrations, the spray RTE for band $n$ becomes:
\be
\bs \cdot \nabla I_{n}(\bx,\bs) =
-\left[\kappa_{{\rm p},n}(\bx) + \overline{\sigma}_{{\rm p},n}(\bx)\right] I_n(\bx,\bs)
+\kappa_{{\rm p},n}(\bx) \; I_{{\rm b,p},n}(\bx)
+\frac{\overline{\sigma}_{{\rm p},n}(\bx)}{4\pi}U_n(\bx)
\ee
where the source function is based on the average particle temperature within a cell.


\subsection{Forward Fraction of Scattering}
\label{forward_fraction}

The effective scattering coefficient in Eq.~(\ref{simple_particle_RTE}) is defined:
\be
\label{effective_sigma_d2}
\overline{\sigma}_{\rm p}(\bx,\la) = \frac{4\pi}{4\pi-\delta\Omega^l} \Big(1-\chi_{\rm f}(\bx,\la)\Big) \; \sigma_{\rm p}(\bx,\la)
\ee
where $\chi_{\rm f}$ is the fraction of the total intensity originally within the solid angle $\delta\Omega^l$ that is scattered into the same angle, $\delta\Omega^l$. The computation of $\chi_{\rm f}$ has been derived in Ref.~\cite{Yang:3}. It can be shown that here $\chi_{\rm f}$ becomes:
\be
\chi_{\rm f}(r,\la) = \frac{1}{\delta\Omega^l}
\int_{\mu^l}^1\int_{\mu^l}^1\int_{\mu_{\rm p,0}}^{\mu_{\rm p,\pi}}
\frac{P_0(\theta_{\rm p},\la)}
{\left[(1-\mu^2)(1-\mu'^2)-(\mu_{\rm p}-\mu\mu')^2\right]^{1/2}}
\; \d\mu_{\rm p} \, \d\mu \, \d\mu'
\ee
where $\mu_{\rm p}$ is a cosine of the scattering angle $\theta_{\rm p}$ and $P_0(\theta_{\rm p},\la)$ is a single droplet scattering phase function
\be
P_0(\theta_{\rm p},\la) =
\frac{\la^2\left(|S_1(\theta_{\rm p})|^2+|S_2(\theta_{\rm p})|^2\right)}{2 \, Q_{\rm s}(r,\la) \pi r^2}
\ee
$S_1(\theta_{\rm p})$ and $S_2(\theta_{\rm p})$ are the scattering amplitudes, given by Mie-theory. The integration limit, $\mu^l$, is a cosine of the polar angle defining the boundary of the symmetric control angle, $\delta\Omega^l$
\be
\mu^l = \cos(\theta^l) = 1 - \frac{2}{N_\Omega}
\ee
The limits of the innermost integral are
\be
\mu_{\rm p,0}   = \mu\mu' + \sqrt{1-\mu^2}\sqrt{1-\mu'^2}  \quad ; \quad
\mu_{\rm p,\pi} = \mu\mu' - \sqrt{1-\mu^2}\sqrt{1-\mu'^2}
\ee


\subsection{Solution Procedure}

The absorption and scattering coefficients, $\kappa_{\rm p}$ and $\overline{\sigma}_{\rm p}$ are not repeatedly calculated during the simulation. Instead, they are tabulated at the beginning of the simulation for each band and a range of different Sauter mean diameters, $r_{32}$. The averaged quantities, now functions $r_{32}$ only, are stored in one-dimensional arrays. During the simulation, the local properties are found by table look-up.

In the band integration of $\kappa_{\rm p}$ and $\overline{\sigma}_{\rm p}$, a constant ``radiation'' temperature, $T_{\rm rad}$, is used to provide the wavelength weighting (Planck function).  $T_{\rm rad}$ should be selected to represent a typical radiating flame temperature. A value of 1173~K is used by default.

The absorption and scattering efficiencies, $Q_{\rm a}$ and $Q_{\rm s}$, and the scattering phase function $P_0(\theta_{\rm p},\la)$, are calculated using the MieV code developed by Wiscombe~\cite{Wiscombe}. The spectral properties of the particles can be specified in terms of a wavelength-dependent complex refractive index. Pre-compiled data are included for water and a generic hydrocarbon fuel based on diesel fuel and heptane. For water, the values of the imaginary part of the complex refractive index (related to the absorption coefficient) are taken from Ref.~\cite{Hale:1}, and a value of 1.33 is used for the real part. For fuel, the droplet spectral properties are taken from Ref.~\cite{Dombrovsky:1}, which includes measurements of the refractive index of a diesel fuel and a comparison of the values with those of heptane. The diesel properties are used for the real part of the refractive index. For the complex part, the heptane properties are used to avoid the uncertainty associated with different types of diesel fuels. Usually, the radiative properties of the particle cloud are more sensitive to the particle size and concentration than to values of the refractive index.



\subsection{Heat Absorbed by Droplets and Small Particles}

The droplet contribution to the radiative loss term is
\be \dot{q}'''_{\rm r} \equiv -\nabla\!\cdot \dbq_{\rm r}''(\bx)(\mbox{droplets}) =
    \kappa_{\rm p}(\bx) \, \left[ U(\bx) - 4 \pi \, I_{{\rm b,p}}(\bx) \right]
\ee
For each individual droplet, the radiative heating/cooling power is computed as
\be
\dq_{\rm r} = \frac{m_{\rm p}}{\rho_{\rm p}(\bx)}
    \kappa_{\rm p}(\bx) \, \left[ U(\bx) - 4\pi \, I_{{\rm b,p}}(\bx) \right]
\ee
where $m_{\rm p}$ is the mass of the droplet and $\rho_{\rm p}(\bx)$ is the total density of droplets in the cell.



\section{Absorption by Large, Non-Scattering Particles}

FDS uses Lagrangian particles to represent objects that are not resolvable on the numerical grid, but are still relatively large compared to liquid droplets. For example, vegetation is modeled as a collection of large, opaque, non-scattering particles. For more discussion of the radiation absorption by large particles, see Section~\ref{rad_part_absorb}.








